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Abstract 

Platyceps rhodorachis (JAN, 1863) and P. karelini (BRANDT, 1838) are both members of the 
P. rhodorachis complex species, which is widely distributed in Iran and includes many local 
populations in the country. These two species are molecularly and to some extent morphologically 
valid. However, hybrids between P. kareliniand P. rhodorachis have been described, but so far their 
ecological differentiation have not been evaluated. In this study, the ecological niche models was 
predicted for these two members of the P. rhodorachis complex using bioclimatic layers and 
geographical coordinates. Possible habitat models show the distribution density of these two species 


in the southern (including some islands in the Persian Gulf) regions, and some areas in northeastern 
Iran. The results of niche similarity tests (identity and niche overlap tests) based on the criteria of 
environmental species, in order to assess the degree of species differentiation, indicate the degree of 
differentiation between these two sister species and raises the possibility of a hybrid belt in southern 
Iran. 
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INTRODUCTION 

The issue of species delimitation has become a contentious issue over the last half century 
(De Queiroz, 2007). Biologists believe that the delimitation of species depends on the identification of 
biological diversity at the species level (Carstens et al., 2013). Numerous operational techniques 
based on molecular or phenotypic data have been designed to aid the understanding of biological 
species, however in some cases post- or pre-zygotic barriers are broken (Mayr 1978; Tobias et al., 
2010). Because selection acts on phenotypes and any mutualistic units, the concept of ecological 
species is proposed as a solution to describe species boundaries (Van Valen, 1976). Interest in 
describing, understanding, and predicting geographic and environmental distributions of species is 
very old (Grinnell, 1917; Wallace, 2016). According to Dobson (1995), "Biodiversity is the diversity 
and variability of living organisms and ecological complexes and the sum of the different types of 
organisms that inherit an area". Understanding how climate change affects the distribution of 
different animal species has always been one of the main challenges for conservation biologists (Root 
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& Schneider, 2006). A set of techniques commonly called species distribution model (SDM) (Guisan & 
Zimmermann, 2000; Guisan & Thuiller, 2005), and the related Ecological Niche Modeling has 
experienced significant growth over the past decade and serves as one of the fundamental strategies 
for determining potential species habitats (Peterson, 2006). The fundamental importance of 
identifying species distribution areas as a biogeographical and ecological concept and the myriad 
applications that can be used through these methods explain the growing interest in this field 
(Peterson et al., 2011). Bioclimatic envelope models is one of the neutral terms that has been 
developed to use these models (Aratijo & Peterson, 2012). 

The concept of ecological species, as mentioned before, has a lot to do with the climate and 
ecosystem of the species distribution region, so the study of species biodiversity in Iran is very 
important due to the existence of different formations in geological, climatic and soil arrangement 
(Agapow et al., 2004; Jowkar et al., 2016). Habitat loss is the greatest threat to reptile and amphibian 
populations, and climate change due to human activities is exacerbating it (Fisher et al., 1968; Devos 
& Wiles, 1973; Ashrafzadeh et al., 2019). Detection of the protective status of reptiles is faced with a 
lack of information for various reasons. Modern computer processing and modeling methods with the 
integration of mechanical models and distribution data of the target species make it possible to gain a 
better understanding of species limiting factors(Root & Schneider, 2006; Yousefkhani, 2019; Ghelichy 
Salakh et al., 2020). The MaxEnt principle is one of the most popular methods for modeling species 
distribution, which is based on ecosystem data and species presence area. The high power of 
presence prediction, low cost and ease of use have attracted many researchers to use the method 
(Hijmans et al., 2005; Warren & Seifert, 2011). 

Colubrid Snakes make up about two-thirds of the world's snakes (Dessauer, 1967). Racer 
snakes of the genus Platyceps distributed from SW Croatia to Central Asia (Kyrgyzstan), Himalayas 
(probably westernmost Nepal) and parts of North and Northeast Africa (Schaetti et al., 2014; Sinaiko 
et al., 2018). So far, 29 species of this genus have been described, of which six species are distributed 
in Iran (Schaetti et al, 2014; Rajabizadeh, 2018). Among Iranian representatives, P. karelini, 
P.ventromaculatus and P. rhodorachis are defined as complex species (Schatti & Mccarthy, 2004). 
Platyceps ventromaculatus (Gray, 1834) is somewhat isolated from P. rhodorachis complex in terms of 
distribution range, but the separation of the distribution range of the other two species has not been 
considered so far (Yildiz, 2011). All research on P. rhodorachis complex has so far focused on 
taxonomic and phenotypic studies and no study has been done on the extent of distribution and 
habitat suitability (Perry, 1985; Khan, 1997; Schatti & McCarthy, 2004; Schaetti et al., 2014; Sinaiko et 
al., 2018). 

In this study, we examined the current distribution of two species (P .rhodorachis and P. 
karelini) using ecological niche criterion, to examine the degree of ecological differentiation between 
the two sister species. Also, modeling was used to predict the potential distribution of both species in 
Iran and the degree of niche space overlap between them. Finally, it is predicted that despite the 
overlap in the distribution range, the focus area of these two species are different from each other 
and we will discuss the effect of non-living environmental factors (temperature, precipitation, 
altitude and slope) on the extent of this separation. Also, we will determine the potential 
hybridization region by using niche modeling. 


MATERIAL AND METHODS 

More of the occurrence records are related to the field studies of the research team of Imam Hossein 
University and Hakim Sabzevari University between 2003 and 2008, which include 126 points. 
Otherwise, the rest are obtained from the literatures (Schatti et al., 2012; Schaetti et al., 2014). In 
total, 158 records of presence belonging to both species (50 records for P. karelini and 108 records 
for P. rhodorachis) (Appendix 1) were used in this study. The current 19 bioclimatic variables were 
loaded from WorldClim (www.worldclim.org) with an accuracy of 30 seconds and used to model both 
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species (Fick & Hijmans, 2017). The slope layer was created using ArcGIS 10.6.1 from the original 
altitude layer in 30 arc-second resolution and all layers were cut using ArcGIS 10.6.1 for Iran. 
Correlation matrix for bioclimatic variables was obtained using ENMtools v1.3 software and Pearson 
correlation coefficient was obtained using ENMTools V 1.3 software and variables with Pearson 
correlation coefficient higher than 0.75 were removed due to high correlation. (Benesty et al., 2009; 
Warren et al., 2010). Finally, three bioclimatic variables and two geographical variables were selected 
for analysis as follows: BIO6 (Min temperature of coldest month of the year (°C), BIO11 (Mean 
temperature of coldest quarter of the year (°C), BIO18 (Precipitation of warmest quarter of the year 
(mm), altitude and slope. 

The species potential distribution model was implemented using the maximum entropy 
method in Maxent v3.4.4 (Elith et al., 2011). The ENMeval package in R v4.0.4 was used to obtain the 
optimal Maxent settings (Muscarella et al., 2014). Maximum 500 iterations, convergence threshold 
was Set to the lowest (0.01) due to the wide range of species distribution, regularization multiplier 
0.5 and 15 replicates with cross-validation method (Heidari, 2019). The area under receiver 
operating characteristic curve (AUC) were considered as accuracy criteria between 0 and 1, the 
amount of lower than 0.5 represents the random model, and a value closer to 1.0 shows high accuracy 
of the predicted model (Townsend Peterson et al., 2007). 

To measure habitat differentiation between species, niche overlap and identity test were 
measured using ASCII files on ENMTools 1.3 (Nakazato et al., 2010). In these tests, two required 
criteria were measured, namely Schoener’s D (Giannini & Goloboff, 2010), and Hellinger’s-based | 
(Schoener & Gorman, 1968). Schoener’s D calculates the suitable range based on the probability of 
occupied grid cells, whereas Hellinger’s-based I refers to the probability of ecological niche models 
(Warren et al., 2010). Identity tests were conducted for each species based on 100 pseudoreplicates. 
Finally, Schoener’s D and Hellinger’s-based I indices of true calculated niche were compared with the 
null distribution of 100 pseudoreplicates (Warren et al., 2008; Gholipur, 2018). These indices ranged 
between 0 (no overlap) and 1 (complete overlap). 


RESULTS 

Ecological niche modeling results 

Based on the literature and records, the distribution range of both species; P. karelini and P. 
rhodorachis overlaps, and there are reports of hybridization among their different populations 
(TyHues, 2000; Schatti et al., 2012). To evaluate the performance of the Maxent model, the area under 
the curve (AUC) was used. The area under the curve (AUC) is a quantitative indicator to show the 
efficiency and accuracy of the model prediction (Elith et al., 2006). The AUC values for P. karelini and 
P. rhodorachis are 0.812+0.166 (mean + standard deviation) and 0.845+0.255 respectively, which is 
acceptable considering the occurrence records. The potential distribution model for P. karelini shows 
the high distribution suitability of this species throughout Iran, except the Zagros and northwestern 
heights of the country; The main distribution for P. karelini is in the southeastern in Sistan and 
Baluchestan Province and to some extent eastern and northeastern of Iran (Figure 1: A). For P. 
rhodorachis, the Niche distribution model shows a high density in the southern coast strip of Iran, 
Khuzestan and western Golestan province (Figure 1: B); both species are mainly distributed in coastal 
strip in the south. BIO6 (Min Temperature of Coldest Month) in P. karelini and altitude in P. 
rhodorachis are the most important factors determining the forecasting model (Table 1). 

Niche overlap between P. karelini and P. rhodorachis indicated that their niche similarity was 
upper than 0.5 (Hellinger’s-based I = 0.89 and Schoener’s D = 0.65) however, due to the hybridization 
between these two sister species, this amount is a significant difference. Identity test showed the 
expected ecological niche interference between two sister species. The results are visible in the 
diagram in Figure 2. The results of the identity test (i = 0.89, D = 0.65) show the relative separation of 
the distribution for these two species, this separation can also be seen in the range of the obtained 
prediction model (Figure 1). 
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TABLE 1. Permutation importance and Percent contribution of variables used in MaxEnt model for 
two species of P. karelini and P. rhodorachis. 


Percent contribution Permutation importance 

Variables Description of variables P.karelini __-P. rhodorachis P. karelini P. rhodorachis 
BIO6 Min Temperature of Coldest Month of 518 20.00 37.4 7.30 

the year (°C) 
BIO11 Mean TeRIDETAure of Coldest Quarter 19.7 32.4 45.7 321 

of the year (°C) 
BIO18 Precipitation of Warmest Quarter of 0.40 050 0.00 1.70 

the year (mm) 
Altitude Height above sea level (m) 26.5 45.6 16.9 58.8 
Slope Land slope 1.60 1.50 0.00 0.00 
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FIGURE 1. Predicted potential distributions of P. karelini (A) and P. rhodorachis (B) using MaxEnt, 
Habitat suitability is shown on the map using a gradient and the white circles indicate specimens 
collecting sites. 
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FIGURE 2. Results of the identity test. Red arrows refer to the actual niche overlap as calculated by 
ENMTools (D and I). The bars (with two different patterns) are calculated by replicates with identity 
test mode. 
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DISCUSSION 

Review of recent papers suggests that ecological theory is rarely explicitly considered (Austin 2007). 
Species distribution models (SDMs), which are commonly used to obtain hypotheses about the 
distribution or realization of species potential, has grown significantly since 1995, when this field 
came into focus. Current linkages between SDM practice and ecological theory are often weak, 
hindering progress. What is clear is that the use of SDM is essential for applications such as 
bioprotection planning and uncertainty in SDMs arises from data deficiencies (eg, samples of species 
occurrences that are small, biased or lacking absences) (Barry & Elith, 2006). Although the 
conservation status of Colubrid snakes has not been studied due to their widespread distribution and 
the least concern for this group has been predicted, but due to the overlap of the distribution range in 
the most suitable habitats of this group with human activities, this group may considered as 
"endangered" in the not too distant future (Petros Lymberakis et al., 2009). Recently, new studies 
have been performed on the P. rhodorachis complex and this group was revised and new species 
recommended (Perry, 2012; Schaetti et al., 2014). Platyceps rhodorachis (Jan, 1865) has the highest 
range of distribution in this genus and overlaps with most congeners and there are significant cases 
of crossbreeding with them, especially with P. karelini (Khan, 1997; TyHues, 2000; Sinaiko et al., 
2018). Based on the literature and records, the distribution range of both species; P. karelini and P. 
rhodorachis overlaps, and there are reports of hybridization among their different populations 
(TyHues, 2000; Schatti et al., 2012), our results confirm this and the range of distribution obtained is 
very similar to reports and literatures. Platyceps karelini (BRANDT, 1838) has three subspecies, all 
three of which are distributed in Iran. There is no report in the literature on the distribution of this 
species in southern Iran and no comprehensive study of the distribution of this species has been done 
so far in Iran. The available studies are limited to the initial faunistic studies and checklists of snakes 
in this area, which due to the observed variations in this genus, there is a possibility of 
misidentification of the species (Rastegar-Pouyani et al., 2011). Comparison of ecological niches 
among different populations can show the niche differentiation and also detect distinct groups along 
the whole distribution range within the species complex (Cayuela et al., 2009; 
Guisan et al., 2013). The predicted distribution pattern in both species has a high overlap with their 
reported distribution range, however areas of the predicted model are reported for the first time. In 
P. karelini, Ardabil and Hormozgan provinces, and especially Sistan and Baluchestan, have been 
identified as compatible areas for the distribution of this species (Figure 1: A). In P. rhodorachis, 
Ardabil province and southwest of Khuzestan province have been identified as areas with high 
compatibility for distribution that Ardabil province had not been reported before (Figure 1: B). 
Altitude is the most important factor determining the likely distribution of P. rhodorachis, while in P. 
karelini the minimum temperature of the coldest month has the greatest effect on the prediction 
model (Table 1). 

In the coastal areas of Sistan and Baluchestan province, the temperature is always above zero 
and due to the high evaporation volume is always sultry and higher than the temperature of the 
province. The same situation is observed in the coastal areas of Hormozgan province, but due to less 
monsoon winds, its intensity has decreased. These conditions justify the desirability of P. karelini 
distribution in these areas (Figure 1) (Chaichitehrani & Allahdadi, 2018). Platyceps rhodorachis 
distribution prediction model indicates the high importance of altitude factor, according to the 
prediction model, this species has the highest distribution merit at low altitudes and close to sea 
level. The precipitation seasonality in Hormozgan province affects the type of vegetation in the region 
(Abolhasan & Maryam, 2013; Nohegar et al., 2015), and this issue directly affects the distribution of 
this species in this area. Probably the distribution of this species in this area is seasonal and reports 
of hybrids between P. rhodorachis and P. karelini in this area confirms this argument (Schatti et al., 
2012). 

The general rule is that the new application of ecological niche modeling greatly facilitates 
species identification, and thus helps to identify additional species diversity and newer species 
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speciation events, it also provides a way to identify the possible hybridization region between the 
two sister species and to explain overlapping distribution patterns. We expect this approach to create 
new opportunities to identify potential distribution areas and protected areas. 
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